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In this work, we explore the use of the semiclassical initial value representation (SC-IVR) method with first- 
principles electronic structure approaches to carry out classical molecular dynamics. The proposed approach can 
extract the vibrational power spectrum of carbon dioxide from a single trajectory providing numerical results that 
agree with experiment and quantum calculations. The computational demands of the method are comparable to 
those of classical single-trajectory calculations, while describing uniquely quantum features such as the zero- 
point energy and Fermi resonances. The method can also be used to identify symmetry properties of given 
vibrational peaks and investigate vibrational couplings by selected classical trajectories. The accuracy of the 
method degrades for the reproduction of anharmonic shifts for high-energy vibrational levels. 
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INTRODUCTION 



Algorithms for the simulation of molecular dynamics be- 
long to the fundamental toolset of modern theoretical chem- 
ical physics. Classical simulation methods are able to study 
systems with up to millions of particles but are unable to de- 
scribe quantum effects such as tunelling and delocalization. 
Exact quantum mechanical methods are restricted to a few 
quantum particles, especially when pre-computed analytical 
potential energy surfaces (PES) are employed. 

First-principles molecular dynamics (FPMD) algorithms 
have been introduced as an alternative to the pre-calculation 
of the PES. FPMD avoids any source of error originated from 
the fitting of the PES. This is particularly true for many de- 
grees of freedom, where the fitting procedure might not repre- 
sent the many-dimensional surface accurately. In this fam- 
ily of methods, the potential and its derivatives are calcu- 
lated on-the-fly as the dynamical simulation progresses and 
are directly obtained from electronic structure calculations. 
In the Born-Oppenheimer molecular dynamics (BOMD) ap- 
proach, the electronic structure calculations for a given sim- 
ulation step are converged based on previous step informa- 
tion. This approach can lead to systematic energy drifts and 
several methods have been proposed to avoid this effect [1]. 
Alternatively, extended Lagrangian molecular dynamics ap- 
proaches (ELMD) [2-5] involve the propagation of nuclear 
and electronic degrees of freedom simultaneously. The elec- 
tronic degrees of freedom are assigned to classical variables 
that are propagated using classical equations of motion and 
these can be expanded in terms of plane waves [2], Gaussian 
functions [4] or real-space grids [5]. Usually ELMD propa- 
gation is computationally more efficient, however questions 
have raised on whether the resultant energy surface remains 
close to the actual Born-Oppenheimer one and about disturb- 
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ing dependencies on the fictitious electronic masses [4, 6]. 

While the evaluation of the potential on-the-fly can be eas- 
ily integrated with classical simulations, the delocalized na- 
ture of quantum mechanical propagation has led to the devel- 
opment of many alternative approaches for the simulation of 
quantum dynamics. For example, the path-integral centroid 
molecular dynamics approach [7] includes quantum nuclear 
effects employing an extended Lagrangian. Alternatively, 
in the variational multi-configuration Gaussian wavepacket 
method (vMCG) [8] the quantum wavepackets are represented 
by fixed-width Gaussian functions for which the potential is 
approximated to be locally harmonic. Other approaches intro- 
duce a mean field approximation and then update the dynam- 
ics in a time-dependent self-consistent fashion [9, 10]. 

Semiclassical molecular dynamics methods [11-20] are 
based on classical trajectories and therefore are amenable 
for carrying out on-the-fly calculation of the potential. The 
benefits of calculating the potential only when needed have 
been suggested by Heller and co-workers [20, 21]. In be- 
tween formally exact quantum methods and classical dynam- 
ics, semi-classical methods include quantum effects approxi- 
mately. Two representative semi-classical approaches are the 
coupled coherent states (CCS) technique [22] and the ab initio 
multiple spawing method (AIMS) algorithm [23]. In the CCS 
approach, several grids of coherent states are classically prop- 
agated and their trajectories can be derived from first principle 
dynamics. In AIMS, the nuclear wavefunction are spawned 
onto a multiple potential surface basis set. This set is made 
of adaptive time-dependent fixed-width Gaussian functions, 
which are generated by classical Newtonian dynamics. 



FIRST-PRINCIPLES SC-IVR 

In this work, we show how the semiclassical initial value 
representation (SC-IVR) [12] method can be coupled tightly 
and naturally, without any mayor change in formulation, with 
first principles electronic structure approaches to carry out 
classical molecular dynamics. We show how the method is 
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able to reproduce approximately quantum effects such as the 
vibrational power spectra using a single, short classical tra- 
jectory using computational resources comparable to those 
employed in first-principles molecular dynamics calculations. 
Calculations employing multiple trajectories can in principle 
be more accurate (and more computational intense as well), 
but here we focus on analyzing the predictive power of single 
trajectory runs. Finally, we describe how different approaches 
can be used in conjunction with this method for studying the 
symmetry of the vibrational states either by arranging the ini- 
tial conditions of the classical trajectory or by employing the 
synnmetry of the coherent state basis. 

In the SC-IVR method, the propagator in F dimension is 
approximated by the phase space integral. 
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where (p (?) , q (f)) are the set of classically-evolved phase 
space coordinates, St is the classical action and Q is a pre- 
exponential factor. In the Heller-Herman-Kluk-Kay [19, 24] 
version of the SC-IVR, the prefactor involves mixed phase 
space derivatives 
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as well as a set of reference states (q | p (f) ,q (f)) = 
Ui ()^/?r)^_/^exp [-7, • {qi - qt (0) /2 + ipi (/) • {qi - qi (/)) /h] 
of fixed width 7. For bound systems, the widths are usually 
chosen to match the widths of the harmonic oscillator approx- 
imation to the wave function at the global minimum and no 
significant dependency has been found under width variation 
[13]. By introducing a 2F x 2F symplectic (monodromy) 
matrix M(t) = ((P(,q«) /"^ (PO)<lo))> one can calculate the 
pre-f actor of Eq. (2) from blocks of F x F size and monitor 
the accuracy of the classical approximate propagation by the 
deviation of its determinant from unity. Wang et al. suggested 
calculating the determinant of the positive-definite matrix 
M-'^M instead [25] and we monitored the same quantity 
for this work. The spectral density is obtained as a Fourier 
transform of the surviving probability [19]. The SC-IVR 
expression of the probability of survival for a phase-space 
reference state \x) = \pN,qN) is 
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The phase-space integral of Eq. (3) is usually computed 
using Monte Carlo methods. If the simulation time is long 
enough, the phase space average can be well approximated by 
a time average integral. This idea has been suggested and im- 



plemented by Kaledin and Miller [26] to obtain the TA (Time 
Averaging [27]) SC-IVR approximation for the spectral den- 
sity. 
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where (p(fi) ,q(fi)) and (p (f2) ,q (fi)) are variables that 
evolve from the same initial conditions but to different times, 
and T is the total simulation time. The advantage of this 
approach is that the additional time integral can in principle 
replace the need for phase-space averaging in the large-time 



limit of a single trajectory. Calculations of the vibrational 
spectra of systems such as the water molecule have proved to 
be very accurate using the TA-SC-IVR approach and its inex- 
pensive single-trajectory variant showed significant improve- 
ments over the simple harmonic approximation for excited vi- 
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brational levels [26]. In order to make Eq. (4) less compu- 
tationally demanding, one can employ the separable approxi- 
mation [26], where the pre-factor of Eq. (4) is approximated 

as a phase, (p (fi ) , q (f i )) = Exp [; (f2) - (f i )) /S] , and 
<^ (f ) =phase[Q(p(0),q(0))]. Using this approximation, Eq. 
(4) becomes 
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leading to a simplification of the double-time integration to 
a single time integral. The resulting integral is positive defi- 
nite, making more amenable for Monte Carlo integration. Our 
numerical tests show that the results of carrying out this ap- 
proximation are essentially identical to the double time inte- 
gral approach when using a single trajectory. In this paper re- 
sults will be reported by use of this last approximation, since it 
is computationally cheaper and numerically more stable than 
Eq. (4). 

For this work, we compute the potential energy surface at 
each nuclear configuration directly from the Kohn-Sham or- 
bitals expanded on a non-orthogonal Gaussian basis. Gradi- 
ents and Hessians at each nuclear configuration are obtained 
analytically from electronic orbitals. The evaluation of the 
potential represents most of the computational effort of our 
approach, which is roughly few hours of computer time using 
standard desktop machines for a 1 cm^' spectrum resolution. 
The nuclear equations of motion are 
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where C is the rectangular matrix of the lowest occupied or- 
bitals and the classical propagation is performed according to 
the velocity- Verlet algorithm, as implemented in the Q-Chem 
package [28]. At each time step, the potential, nuclear gradi- 
ent and Hessian are used to calculate the action, pre-factor and 
coherent state overlaps necessary for the TA-SC-IVR method 
(Eqs. 4 and 5). A schematic representation of an implemen- 
tation of the algorithm for a multithreaded machine is shown 
in Fig. (1). At each time step, results are accumulated for 
time-average integration. The results presented on this work 
were carried out on a single thread. For each classical trajec- 
tory, the procedure is repeated and the final integration gives 
the spectrum intensity / {E) for a given parametric value of 
E. The same procedure is repeated for next E + h.E, where in 
our calculation A£ = 1cm"'. As previously mentioned, the 
trajectory is monitored by calculating at each time step the 
deviation of the determinant of the monodromy matrix from 
unity. The difference in the determinants was always smaller 
than 10~^ during the course of the calculations. A time step of 
10 a.u. has been always found to satisfy the strict monodromy 
matrix restrictions even for the lightest atoms. 




Figure 1 : First-principles SC-IVR algorithm: At each time step elec- 
tronic wavefunction are saved to calculated nuclear Hessian. Nu- 
clear positions, gradients and Hessian are accumulated for the spec- 
tral time-average integral. 



The calculation of the full dimensional vibrational power 
spectrum of the CO2 molecule is a challenging test for FP-SC- 
IVR method: A successful method should reproduce spectral 
features such as degenerate bending modes, strong intermodal 
coupUngs and Fermi resonances. To evaluate the FP-SC-IVR 
method, we compare vibrational spectrum of CO2 molecule 
from FP-SC-IVR method to numerically-exact discrete vari- 
able representation (DVR) eigenvalue calculations on a po- 
tential fitted to a set of first-principles points obtained at the 
same level of theory. The next section describes the details 
of the potential fitting and DVR calculation. Following, we 
continue on the discussion of the FP-SC-IVR method. 



POTENTIAL FITTING AND GRID CALCULATIONS 

The CO2 molecule is a linear molecule with four vibrational 
normal modes: a symmetric stretching mode (Vi), degenerate 
bending modes (V2 and V2) , and an antisymmetric stretching 
mode (V3). A 3c/ potential energy grid in internal coordinates 
is calculated using the B3LYP density functional [29] with 
the cc-pVDZ basis set [30]. The grid points are then fitted to 
a potential energy surface [31] represented by a fourth-order 
Morse-cosine expansion 



X {cose - cosOeY (1 - e-«2('-2-'-''))'' (7) 



where the parameter = 2.206119 a.u. and de = 180 specify 
the equilibrium coordinates of the CO2 molecule. The Morse 
parameters ai = 02 = 1.2489 a.u. were determined so as to 
minimize the standard deviation of the differences of the fit- 
ted potential from the ab initio result using the Levenberg- 
Marquardt non-linear least square algorithm [32] . Instead, rg 
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Table I: Expansion coefficients Kjji^ for the CO2 B3LYP/cc-pVDZ 
fitted potential energy surface in atto Joule units. 



was obtained by geometry optimization within the Q-Chem 
ab initio package [28]. 

The 35 Kijk coefficients were subject to the non-linear least 
square fitting procedure to the DFT energies. Since these co- 
efficients must be the same once rj and r2 are swapped, 13 
linear constraints of the type Kjji^ = K^ji were imposed during 
the fitting procedure. Additionally, to ensure that the equi- 
librium geometry was fitted to the predetermined equilibrium 
parametric distance, the coefficients Kiqq and Afloi were con- 
strained to be zero. Consequently, we employed a total num- 
ber of 14 fitting constraints {Kqqq term is always constant). A 
total of 2500 ab initio grid points were chosen for the fitting 
process. These grid points range from 1.42 a.u. to 7.09 a.u. 
for ri and r^, and from 1 13.6 to 180 for the angle variable. The 
calculated expansion coefficients Kij^ are reported in Tab. (I). 

As far as the numerically exact eigenvalues calculations 
is concerned, we used an exact DVR (Discrete Variable 
Representation) matrix diagonalization procedure. The CO2 
molecule was described for grid calculations in internal co- 
ordinates, while on-the-fly classical trajectories and the semi- 
classical calculations described previously were performed in 
Cartesian coordinates. No significant contamination between 
the rotational (set to zero kinetic energy) and vibrational mo- 
tion was found within the simulation time. To this end, the de- 
viation from simplecticity of the monodromin matrix in the vi- 
brational sub-space were never more than 10~^ as previously 
mentioned. 

The coordinates ri and r2 are CO distances, and is the 
angle between the CO bonds. In these coordinates the kinetic 
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The carbon mass were taken to be mc = 12.0 a.m.u., while the 
oxygen mass mo = 15 .9949 a.m.u. and the reduced mass is as 
usual 1 / jXco = 1 /'"c + 1 / 'Wo- 

As previosuly mentioned, in order to calculate exact eigen- 
values, a sine-DVR basis for the coordinates ri and r2 and a 
Legendre-DVR basis for 6 has been used [33]. For each de- 
gree of freedom 50 DVR functions were used and eigenvalues 
were converged to at least 10~^cm~'. The sine-DVR ranged 
from 1.51 a.u. to 3.78 a.u. and the magnetic quantum number 
m of the Legendre-DVR was zero. 

Because of the restriction of total angular momentum / = 0, 
we couldn't observe all degenerate bending excitations. How- 
ever, ZPE and several vibrational energy levels were obtained 
and compared with that ones coming from a single on-the-fly 
semiclassical trajectory. 



FIRST-PRINCIPLES SC-IVR CALCULATIONS 

The full power spectrum obtained using Eq. (4) after 
3000 BOMD steps of 10 a.u. each is shown on the bot- 
tom of Fig. 2. For longer simulations, the monodromy ma- 
trix symplectic properties as well as the resolution of the 
spectrum started to deteriorate. The calculated vibrational 
zero-point energy (ZPE) value was 2518 cm"' versus the 
exact value of 2514.27 cm~^ and both are in good agree- 
ment with the experimental value of 2508 cm"'. In con- 
trast, harmonic normal-mode analysis (whose frequencies are 
656.62, 1 363 .46, 2423 .47 wavenumbers) predicts a frequency 
of 2550.08 cm"'. Thus, the TA-SC-IVR method successfully 
reproduces the ZPE anharmonic effects with the use of a sin- 
gle classical trajectory. Some representative frequencies of 
the power spectrum are presented in Table II. The ZPE was 
shifted to zero for comparison with reported classical ELMD 
simulations on the same system that cannot reproduce the ZPE 
or higher vibrational states [34, 35] but only single modes fre- 
quencies. For these studies of Refs. [34, 35], the vibrational 
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Figure 2: CO2 Vibrational Power Spectrum: Initial kinetic energy on: 
(a) all modes; (b) symmetric mode; (c) one bending and symmetric 
modes; (d) bending and asymmetric modes. 

data were obtained from the Fourier transform of correlation 
functions of classical trajectories in plane-wave DFT calcu- 
lations. The ELMD approach predicts the following funda- 
mental frequencies 648, 1368, 1428 and 2353 for Ref. [34] 
and 663, 1379, 1456 and 2355 for Ref. [35]. These classical 
results are similar but limited to a normal mode analysis. 

Table 11 compares our TA-SC-IVR results with the exact 
ones and to those obtained by Filho [36] with the same density 
functional and a basis set of comparable quality (6-31H-G*) 
[37], using a perturbative approximation of the eigenvalue ex- 
pansion. One can see how a different basis set results a sig- 
nificant deviation of vibrational levels spacing, once the com- 
parison is performed in units of wavenumbers. 

A major difficult on the CO2 power spectrum simulations 
is the calculation of the Fermi resonance splittings. These are 
the result of anharmonic couplings, and they represent a strin- 
gent test for a semi-classical method that reUes on a single 
short trajectory. The Fermi resonances occur when an acci- 
dental degeneracy between two excited vibrational levels of 
the same symmetry exists and it results in a repulsion between 
the corresponding energy levels. The sources of these reso- 
nances are purely anharmonic and are only present in poly- 
atomic potentials. For the CO2 molecule, the unperturbed fre- 
quencies for the symmetric stretching are roughly equal to the 
first bending overtone (Vj = 2V2). For these modes, the wave- 
functions are transformed as the irreducible representation of 
D„h, i.e. Vi (10*^0) as E+, at the experimental frequency of 
1388 cm" \ and v|(02°0) as E+ + A^, at an experimental fre- 
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"Experimental frequencies in cm" from Ref. [41] 

*First number is the symmetric stretch quantum, second are the degenerate 
bendings, and third one is the asymmetric stretch. The exponent of the second 
number is the /, degeneracy index. 

■^Vibrational levels according to a normal modes harmonic model 

''Using the Separable approximation of Eq.(5) 

Table II: Some of the calculated vibrational energy eigenvalues. All 
data are in wavenumbers. Fermi Resonances group of frequencies 
are indicated by the same superscript symbols. Uncertain peaks are 
marked with (+). The first column represents the experimental vi- 
brational frequencies associated with the modes listed on the second 
column. The third column shows the harmonic DFT results. In the 
fourth and fifth columns, we show our FP-SCIVR and exact numeri- 
cal DVR calculations in the B3LYP/cc-PVDZ model chemistry used 
for the FP-SCIVR calculations. The fifth column shows perturbative 
DFT calculations carried out using a similar functional and basis set. 

quency of 1285 cm" ' . Another Fermi doublet results from the 
addition of a quantum of bending mode to the previous Fermi 
doublet to yield the following states: ViV2(ll'0) , at an ex- 
perimental frequency of 2077 cm"' and the V2(03'0) state, 
at an experimental frequency of 1932 cm" ^ Higher-energy 
Fermi resonances are indicated in Table II by using the same 
superscript symbols. The first Fermi terms are located at 1313 
and 1363 in a harmonic approximation and corrected to 1288 
and 1381 wavenumbers for FP-TA-SC-IVR. Thus, the origi- 
nal levels have been repelled by Fermi couplings. One mode 
is located at a higher frequency than the harmonic prediction, 
while the other is at a lower frequency. The latter effect could 
be explained also by simple anharmonicity, but the former 
is evidence of the ability of the single trajectory FP-TA-SC- 
IVR method even when the separable approximation is used 
to capture Fermi resonance effects partially. The same rea- 
soning can explain the second Fermi doublet located at 1932 
and 2106 for FP-TA-SC-IVR, while the harmonic estimate at 
1970 and 2020 wavenumbers. 

With the FP-TA-SC-IVR method, one can also identify the 
couplings between vibrational modes and the appearance of 
Fermi resonance splittings by carrying out simulations with 
different initial conditions. This can be achieved by selec- 
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tively setting the initial velocity of some vibrational modes to 
zero. The anharmonic coupling between levels leads to a con- 
sistent reproduction of the ZPE peak in the spectrum for all 
simulations. However the excited vibrational peaks related to 
the modes with zero initial kinetic energy show a very small 
signal in the power spectrum. Vibrational energy redistribu- 
tion processes can be studied as well, by carrying out simula- 
tions at different timescales. In Fig. 2, we show the resulting 
power spectra for different initial conditions. If the initial state 
contains only purely symmetric motion, the lowest Fermi res- 
onance peaks in Fig. 2(b) are absent as well as for a bending 
(without symmetric stretching) motion in Fig. 2(d). These 
results and the intensity of their peaks respect to that ones lo- 
cated at the same frequencies in Fig. 2(a) suggest that the 
Fermi resonance is indeed originated from the coupUng be- 
tween bending and the symmetric modes. One can reach the 
same conclusions by inspecting the lower Fermi doublet peaks 
intensity: by adding a bending mode (from Fig. 2(b) to Fig. 
2(c)) and a second one (from Fig. 2(c) to Fig. 2(a)) the inten- 
sity of both peaks is gradually raised. This is called "intensity 
borrowing" and it arises from the strong mixing of the zero 
order states. These observations reinstate that "repulsion and 
mixing are the hallmarks of Fermi resonances" [38]. Also, for 
a distinct set of initial conditions, an additional peak at 5500 
cm^ ' related to the asymmetric stretch was observed. Using 
the proposed approach, one can carefully detect the character- 
istics of each peak even for complicated power spectra. 

An attractive method for obtaining the symmetry prop- 
erties of the eigenstates involves arranging the initial ba- 
sis vectors [26, 39]. The basis for this method is the di- 
rect product of coherent states \x) = I\t=i\pf\^f^)^'' ■ 
These states can be chosen to have an initial symmetry by 
employing linear combinations of the form \pf\qf''Y'' = 
{\pf\q'h+£k\-p,-qT')) /\/2. The A;-th mode can be 
made synometric (e*: = 1), antisymmetric (e;t = — 1) or have no 
symmetry restrictions (e^. = 0). In order to assign the proper 
symmetry to each peak on Fig. 3 , the reduced Djh symme- 
try group was adopted. AH irreducible representations were 
reproduced and peaks were grouped by symmetry as reported 
in Fig. 3. Note that (d) and (e) plots are identical since they 
only differ trivially by swapping coefficients between the de- 
generate bending modes in the original Dooh symmetry group. 

Finally we have investigated the stability of the propagator 
versus variations of the coherent states gaussian width param- 
eters jj. Previous calculations [24] have shown that there is no 
significant depedency on energy and norm conservation for 
the semiclassical propagator if suitable values of ji are cho- 
sen. For power spectra calculation we have chosen to look 
at vibrational levels variations under different values of co- 
eherent states width. Since a single trajectory was used in the 
FP-TA-SC-IVR approach, no Monte Carlo integration is per- 
formed in phase space coordinates and the changes of ji are 
confined to the coherent states overlap and to the prefactor in 
Eq. (2). As reported in Fig. 4 and checked on a finer scale, no 
significant variation was observed beyond 1 cm" ^ . These find- 
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Figure 3: CO2 Vibrational Power Spectrum (Separable approxi- 
mation): Different basis set symmetries for Vj (symmetric stretch- 
ing mode), V2 and V2 (bending modes) and V3 (asymmetric mode) 
and the corresponding Dih irreducible representation: (a) all es are 
zero: (b) (5i„): e(vi) = 0,£(v2) = l,e(v2) = 0,£(v3) = -1; (c) 
{Ag): e(vi) = 1,£(V2) = 0,e(V2) =0,e(v3) = 1; (d) (S2„):£(vi) = 
0,£(V2) = -l,e(v2) =0,£(v3) = 1, (e) (S3J £(v,) = 0,e(v2) = 
0,e(v2) = — l,e(v3) = 1. B2u and63„ representations are degener- 
ated in the D„h subspace as shown. 



ings are in agreements with previous calculations on the same 
propagator [24]. Interestingly, a different distribution in peaks 
intensity were found in each panel. Since the peaks magnitude 
is proportional to the overlap between the reference state and 
the actual eigenfunction, the anharmonic choice (7; = co,/2) is 
a more suitable solution as clearly showed on panel (c) of Fig. 
4. 



CONCLUSIONS 

In conclusion, we have shown that SC-IVR can be imple- 
mented easily and efficiently using first principles molecular 
dynamics. With the modest computational cost of a single 
classical trajectory, the vibrational density of states of the CO2 
molecule was calculated. On Fig. 5 we report a graphical 
comparison between the harmonic and the FP-TA-SC-IVR ap- 
proximations, versus the exact vibrational value for the Fermi 
resonance multiplets. One can notice how the single trajec- 
tory FP-TA-SC-IVR goes far beyond the harmonic approxi- 
mation by removing the harmonic degenerancy and including 
part of anharmonicity. Fermi splittings are well mimiced not 
only for the first doublet, but also for the higher ones. The 
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Figure 4: Gaussian width variations and related power spectra: a) 
){ = CO,-; h)Yi = 2cOi\ c)yi = (Oi/2, where to, are the i — esime normal 
mode frequency. The FP-SCIVR power spectra are fairly insensitive 
to variations in the value of the coherent state width. 
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Figure 5: Fermi Resonance states vibrational energy level: (a) in har- 
monic approximation; (b) single FP-SC-IVR trajectory calculation; 
(c) exact grid calculation on spUned potential. 



numerically exact DVR vibrational energy levels constrained 
by / = are represented on the last column. The FP-TA- 
SC-IVR values are similar to the DVR results, when compar- 
ison is possible. However, a closer look at Table (II) shows 
how these single trajectory FP-TA-SC-IVR calculations can 
include only part of the anharmonicity and that their preci- 
sion gets worse for higher vibrational levels. In particular, the 
spacing of the higher-energy states is harmonic-like and thisis 
the mayor limitation of using a single classical trajectory. 

These and previous calculations on model potentials [26] 
has shown how the single trajectory TA-SC-IVR gives rea- 
sonable results and performs better for higher frequencies 
modes. The computational cost of the method is essentially 
the same as classical propagation, and therefore, if broadly 
implemented in electronic structure codes, it can provide a 
description of quantum effects at a comparable computational 
cost to that of classical approaches. Possible applications of 
this method or related ones are the study of excited electronic 
states and Franck-Condon transitions, such as vibrational ab- 
sorption spectra [42]. Although this single trajectory approach 
may be a practical tool for the simulation of more complex 
systems, the use of more trajectories is probably required to 
remove any harmonic"ghost states". We are currently explor- 
ing the use of a small number of a set of systematically de- 
termined trajectories for further improvement of the results. 
If the number of required trajectories grows as a low polyno- 
mial of the system size, semi-classical methods could be com- 
petitive with currently-employed numerical approximations to 



obtain anharmonic vibrational effects. Finally, we expect that 
the representation of the potential energy in terms of normal 
coordinates will become less suitable when large amplitude 
motions or non adiabatic effects come into play. 



Acknowledgement 

One of the authors (M.C.) feels deeply in debt with Prof. W. 
H. Miller for the many lessons learned from him. The authors 
thank Dr. A. Kaledin, Prof. E. J. Heller for useful discussions 
and revision of the manuscript. A. A.-G., S. S. and S. A. thank 
the Faculty of Arts and Sciences of Harvard University for 
financial support and S. S. thanks the Samsung Scholarship for 
financial support. M. C. and G. F. T. thanks the University of 
Milan for fundings and CILEA (Consorzio Interuniversitario 
Lombardo per L'Elaborazione Automatica) for computational 
time allocation. A. A.-G., S. S. and S.A thank FAS Research 
Computing for cluster computing support. 



' Electronic address: michele . ceottoQunimi . it , aspuruQ 
chemistry . harvard . edu 
[1] J. M. Herbert and M. Head-Gordon, Phys. Chem. Chem. Phys., 
2005, 7, 3269. 

[2] R. Car and M. Parrinello, Phys. Rev. Lett, 1985, 55, 2471. 
[3] H. B. Schlegel, J. M. Millam, S. S. Iyengar, G. A. Voth, A. D. 
Daniels, G. E. Scuseria, and M. J. Frisch, / Chem. Phys., 2001, 



8 



114, 9758. 

[4] J. M. Herbert and M. Head-Gordon, / Chem. Phys., 2004, 121, 
11542. 

[5] Y. Liu, D. Yame, M. E. Tuckerman, Phys. Rev. B, 2003, 68, 
125110. 

[6] P. Tangney, J. Chem. Phys., 2006, 124, 044111. 
[7] M. Pavese, D. R. Berard, G. A. Voth, Chem. Phys. Lett, 1999, 
300, 93. 

[8] G. A. Worth, M. A. Robb and I. Burghardt, Faraday Discuss., 
2004, 127, 307. 

[9] S. Iyengar and J. Jakowski, J. Chem. Phys., 2005, 122, 114105. 
[10] O. Knospe and P Jungwirth, Chem. Phys. Lett., 2000, 317, 529. 
[11] W. H. Miller, Adv. Chem. Phys., 1974, 25, 69; W. H. Miller, 

Faraday Discuss., 1998, 110, 1. 
[12] W. H. Miller, / Chem. Phys., 1970, 53, 3578; ibidem, 1970, 

53, 1949; W. H. Miller, J. Phys. Chem. A, 2001, 105, 2942; M. 

Thoss and H. Wang, AnnM. Rev. Phys. Chem., 2004, 55, 299; K. 

G. Kay, Annu. Rev. Phys. Chem., 2005, 56, 255. 
[13] H. Wang, X. Sun, and W. H. Miller, J. Chem. Phys., 1988, 108, 

9726; X. Sun and W. H. Miller, J. Chem. Phys., 1999, 110, 

6635; M. Thoss, H. Wang and W. H. Miller, J. Chem. Phys., 

2001, 114, 9220; T. Yamamoto, H. Wang, and W. H. Miller, J. 

Chem. Phys., 2002, 116, 7335; T. Yamamoto, W. H. Miller, J. 

Chem. Phys., 2003, 118, 2135. 
[14] J. Ankerhold, M. Saltzer, and E. Pollak, J. Chem. Phys., 2002, 

116, 5925; S. Zhang and E. Pollak, Phys. Rev Lett, 2003, 91, 

190201; S. S. Zhang and E. Pollak, J. Chem. Phys., 2004, 121, 

3384. 

[15] A. R. Walton, D. E. Manolopoulos, Mol. Phys., 1996, 87, 961; 

A. R: Walton, D. E. Manolopoulos, Chem. Phys. Lett., 1995, 

244, 448; M. L. Brewer, J. S. Hulme, D. E. Manolopoulos, J. 

Chem. Phys., 1997, 106, 4832. 
[16] S. Bonella, D. Montemayor , and D. F. Coker, Proc. Natl. Am. 

Soc, 2005, 102, 6715; S. Bonella and D. F. Coker, J. Chem. 

Phys., 2003, 118, 4370. 
[17] Y. Wu , M. Herman, V. S. Batista, J. Chem. Phys., 2005, 122, 

114114; Y. Wu and V. S. Batista, J. Chem. Phys., 2003, 118, 

6720. 

[18] F. Grossmann, Comments At. Mol. Phys., 1999, 34, 243. 

[19] E. J. Heller, J. Chem. Phys., 1975, 62, 1544; E. J. Heller, J. 

Chem. Phys., 1981, 75, 2923. 
[20] E. J. Heller, Acc. Chem. Res., 1981, 14, 368; E. J. Heller, Acc. 

Chem. Res., 2006, 39, 127. 
[21] T. Van Voorhis, and E. J. Heller, J. Chem. Phys. , 2003, 119, 

12153 

[22] D. V. ShalashiUn and M. S. Child, Chem. Phys., 2004, 304, 



103; D. V. ShalashiUn and M. S. Child, J. Chem Phys., 2001, 
115, 5367. 

[23] M. Ben-Nun, T. J. Martinez, Adv Chem. Phys., 2002, 121, 439. 

[24] M. R Herman and E. Kluk, Chem. Phys., 1984, 91, 27; K. G. 
Kay, /. Chem. Phys., 1994, 100, 4377; K. G. Kay, J. Chem. 
Phys., 1994, 100,4432. 

[25] H. Wang, D. E. Manolopoulos, and W. H. Miller, J. Chem. 
Phys., 2001, 115, 6317. 

[26] A. L. Kaledin and W. H. Miller, J. Chem. Phys., 2003, 118, 
7174; M. Ceotto, PhD Dissertation, University of California, 
Berkeley (2005); A. L. Kaledin and W. H. Miller, / Chem. 
Phys., 2003, 119, 3078. 

[27] Y. Elran and K. G. Kay, J. Chem. Phys., 1999, 110, 3653; ibi- 
dem, 1999, 110, 8912. 

[28] Y. Shao, et al. Phys. Chem. Chem. Phys., 2006, 8, 3172. 

[29] A.D. Becke, / Chem. Phys., 1993, 98, 5648; R J. Stephens, R 
J. Devlin, C. F. Chabalowski, and M. J. Frisch, /. Phys. Chem., 
1994,98, 11623. 

[30] T. Dunning, Jr. J. Chem. Phys., 1989, 90, 1007. 

[31] J. Zuniga, M. Alacid, A. Bastida, F. J. Carvajal, and A. Re- 
quena, / Mol. Spectr, 1999, 195, 137. 

[32] K. Levenberg, Quart. Appl. Math. 2, 164 (1944); D. Mar- 
quardt, Siam J. Appl. Math. 11, 431 (1965); M.I.A. Lourakis, 
Levenberg-Marquardt nonlinear least squares algorithms in 
C/C-l~l-, http://www.ics.forth.gr/~lourakis/levmar/ (2004) 

[33] M. H. Beck and H.-D. Meyer, J. Chem. Phys., 2001, 114, 2036; 
G. A. Worth, M. H. Beck, A. Jackie, and H.-D. Meyer, The 
MCTDH Package, Version 8.2, (2000), University of Heidel- 
berg, Heidelberg, Germany, H.-D. Meyer, Version 8.3 (2002). 
See http://www.pci.uni-heidelberg.de/tc/usr/mctdh/ 

[34] R Gygi, Phys. Rev. B, 1995, 51, 11190. 

[35] J. R. Chelikowsky, X. Jing, K. Wu, and Y. Saad, Phys. Rev. B, 
1994, 53, 12071. 

[36] H. R M. Filho, Spectr Acta Part A, 2002, 58, 2621. 

[37] W. J. Hehre, R. Ditchfield and J. A. Pople, J. Chem. Phys., 1972, 
56, 2257. 

[38] E. J. Heller, E. B. Stechel, and M. J. Davis, J. Chem. Phys., 
1980, 73, 4720. 

[39] X. Sun and W. H. Miller, J. Chem. Phys., 1998, 108, 8870. 

[40] A. M. N. Niklasson, C. J. Tymczak, and M. Challacombe, Phys. 
Rev Utt, 2006, 97, 123001; A. M. N. Niklasson, C. J. Tym- 
czak, and M. Challacombe, /. Chem. Phys., 2007, 126, 144103. 

[41] L. R. Brown and C. B. Farmer, Appl. Opt., 1987, 26, 5154. 

[42] J. Tatchen and E. Pollak, J. Chem. Phys., 2009, 130, 041103. 



